Average Bubble Profile

In [1]:
import numpy as np
import math
import matplotlib
#matplotlib.rcParams.update({'font.size': 12})
from matplotlib import gridspec
import matplotlib.pyplot as plt
from collections import deque
import scipy as scp
import scipy.optimize as sco
import scipy.signal as scs
import scipy.special as ssp
import scipy.integrate as sci
import scipy.interpolate as intp
import statistics as stat
import random
from functools import partial
from operator import eq
from itertools import zip_longest, compress, count, islice, groupby, cycle
from labellines import labelLine, labelLines
import os

np_load_old = np.load
np.load = lambda *a,**k: np_load_old(*a, allow_pickle=True, **k)

Simulation Parameters

In [2]:
nLat = 4096
nSims = 390
minSim = 0
simstep = 1

nu = 2.*10**(-3)
lamb = 6; print('lamb = ', lamb)
m2eff = 4. * nu * (- 1. + lamb**2); print('m2eff = ', m2eff)
lenLat = 2 * 50. / np.sqrt(2. * nu); print('lenLat = ', lenLat)
phi0 = 2. * np.pi / 7.; print('phi0 = ', phi0)

alpha = 8.
nCols = 1
phi_initial = np.pi
dt_phi_initial = 0.
mask = 4*phi_initial

############################################################
nyq = nLat//2+1
spec = nyq
hLat = nLat//2
dx = lenLat/nLat; print(dx)
dk = 2.*np.pi/lenLat
dt = dx/alpha
dtout = dt*alpha
light_cone = dtout/dx
lamb =  6
m2eff =  0.28
lenLat =  1581.1388300841897
phi0 =  0.8975979010256552
0.3860202221885229
In [3]:
titles = [r'$\phi(x)$', r'$\partial_t \phi(x)$', r'$|\nabla \phi(x)|^2$', r'$V(\phi(x))$']
plots_file = '/home/dpirvu/big_plot_file/thin_wall_average_bubble/'
pickle_location = '/home/dpirvu/pickle_location/thin_wall_average_bubble/'
cutout_location = '/gpfs/dpirvu/thin_wall_average_bubble/'
suffix = '_for_phi0{:.4f}'.format(phi0)+'_len{:.4f}'.format(lenLat)+'_lamb{:.4f}'.format(lamb)+'_x'+str(nLat)

def bubbles_file(min, max):
    return cutout_location+'bubbles_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def centred_bubbles_file(min, max):
    return pickle_location+'centred_bubbles_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def bubbles_at_rest_file(sim):
    return cutout_location+'restbubble_sim'+str(sim)+suffix+'.npy'
def average_bubble(min, max):
    return pickle_location+'average_bubble_from_sim'+str(min)+'_up_to_sim'+str(max-1)+suffix+'.npy'
def average_of_N_bubbles(Nbubbles):
    return pickle_location+'average_of_'+str(Nbubbles)+'_bubbles'+suffix+'.npy'
In [4]:
instanton_location = '/home/dpirvu/inst/instantons/dev/thin_wall_instanton_sim.dat' #/np.sqrt(m2eff)/dx
a = np.genfromtxt(instanton_location)
coleman_profile = np.pi+a[:,0]
xoffset = nLat//8
temp = coleman_profile[len(coleman_profile)//2-xoffset:len(coleman_profile)//2+xoffset]
peaks, _ = scs.find_peaks(temp)
fwhm, height, left_ips, right_ips = scs.peak_widths(temp, peaks, rel_height=0.5)
radius_Coleman_bubble = max(fwhm)/2.
filter_size = radius_Coleman_bubble*dx
print('filter_size = ', filter_size)

plt.plot(np.arange(len(temp))-xoffset, temp)
plt.plot(peaks-xoffset, temp[peaks], "x")
[plt.hlines(height[i], left_ips[i]-xoffset, right_ips[i]-xoffset, color="C2") for i in range(len(fwhm))]
plt.xlabel(r'$r_{\rm E}$'); plt.ylabel(r'$\phi$')
filter_size =  67.08180510014577
/cm/shared/apps/conda-environments/python37/lib/python3.7/site-packages/ipykernel_launcher.py:7: PeakPropertyWarning: some peaks have a width of 0
  import sys
Out[4]:
Text(0, 0.5, '$\\phi$')
In [5]:
def V(phi):
    return ( -np.cos(phi) + 0.5 * lamb**2. * np.sin(phi)**2. ) * 4. * nu
def dV(phi):
    return ( np.sin(phi) + 0.5 * lamb**2. * np.sin(2.*phi) ) * 4. * nu

right_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[np.pi, 2*np.pi], method='bounded')
left_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[0, np.pi], method='bounded')
right_left_phi_at_V_max = sco.minimize_scalar(lambda x: -V(x), bounds=[2*np.pi, 3*np.pi], method='bounded')

def F(x):
    return V(x) - V(phi_initial)
phi_upper_bound = sco.fsolve(F, 5.5)[0]
phi_lower_bound = sco.fsolve(F, 0.5)[0]
phi_lower_bound2 = sco.fsolve(F, -0.5)[0]
print(phi_lower_bound2)
phi_upper_lower_bound = sco.fsolve(F, 6.5)[0]

############################################################
wall_tension, err = sci.quad(lambda x: np.sqrt(2*(V(x) - V(2*np.pi))), np.pi, 2*np.pi); print('wall_tension = ', wall_tension)
epsilon = V(np.pi) - V(2*np.pi); print('epsilon = ', epsilon)
R_coleman = wall_tension/epsilon; print('R_coleman = ', R_coleman)
wall_thkn = right_phi_at_V_max.x / np.sqrt( V(right_phi_at_V_max.x) - V(np.pi) ); print('wall_thkn = ', wall_thkn)
############################################################

matplotlib.rcParams.update({'font.size': 12})
fig = plt.figure(figsize=(7, 4))
plt.plot([i for i in np.arange(0, mask, phi_initial/100)], [V(i) for i in np.arange(0, mask, phi_initial/100)])
plt.plot(phi_upper_bound, V(phi_upper_bound), 'ro')
plt.plot(phi_lower_bound, V(phi_lower_bound), 'co')
plt.plot(phi_lower_bound2, V(phi_lower_bound2), 'co')
plt.plot(phi_upper_lower_bound, V(phi_upper_lower_bound), 'yo')

plt.plot(right_left_phi_at_V_max.x, V(right_left_phi_at_V_max.x), 'ko')
[plt.plot(i*phi_initial, V(i*phi_initial), 'bo') for i in range(5)]
[plt.axvline(left_phi_at_V_max.x + i*phi_initial, color='y', ls='--') for i in range(0,3,2)]
[plt.axvline(right_phi_at_V_max.x + i*phi_initial, color='orange', ls='--') for i in range(0,3,2)]
plt.xlabel(r'$\phi(x)$'); plt.ylabel(r'$V(\phi(x))$'); plt.show()
-0.33489615843937864
wall_tension =  1.162408162505219
epsilon =  0.016
R_coleman =  72.65051015657619
wall_thkn =  12.697740180727095

VISUALISING SIMULATION DATA

In [6]:
def add_mask(field_slice, threshold):
    return field_slice * [0 if np.abs(i) > threshold else 1 for i in field_slice]

def smoothen(field_slice, sigma, bool):
    N = len(field_slice)
    if bool:
        pbc = [j-1 if j < N/2+1 else N-j+1 for j in range(1, N+1)]
    else:
        pbc = [j-1 for j in range(1, N+1)]
    window = [np.exp(- 0.5 * (x*dx/sigma)**2) / np.sqrt(2*np.pi) / (sigma/dx) for x in pbc]
    window = window / sum(window)
    spectral_filter = np.fft.fft(window, len(field_slice))

    fft_field_slice = np.fft.fft(field_slice, len(field_slice))
    smooth_fft_field_slice = [spectral_filter[k] * fft_field_slice[k] for k in range(len(fft_field_slice))]
    smooth_field_slice = np.fft.ifft(smooth_fft_field_slice, len(field_slice))
    return np.asarray([i.real for i in smooth_field_slice])

def find_slice_peaks(field_slice, peak_threshold):
    """ Finds x coordinate of peaks in masked field with mask applied at threshold. """
    peak_coord = signal.find_peaks(field_slice, height = peak_threshold)[0].tolist()
    if field_slice[-1] > peak_threshold and field_slice[0] > peak_threshold and field_slice[-1] != field_slice[0]:
        if field_slice[0] > field_slice[-1] and field_slice[0] > field_slice[1]:
            peak_coord.append(0)
        elif field_slice[0] < field_slice[-1] and field_slice[-1] > field_slice[-2]:
            peak_coord.append(len(field_slice)-1) # this minds potential boundary discontinuities
    peak_heights = [field_slice[coord] for coord in peak_coord]
    return peak_coord, field_slice.tolist().index(np.max(peak_heights))

def truncateNum(num, decimal_places):
    StrNum = str(num)
    p = StrNum.find(".") + 1 + decimal_places
    return float(StrNum[0:p])

def time_at_fraction(bubble, frac, limit):
    T, N = len(bubble), len(bubble[0])
    right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
    time_list = [t if (right_phi_x[t] <= N*frac) else 0 for t in range(T)]
    return next((t for t in time_list[::-1] if t != 0), 0)

def time_at_size(bubble, size, limit):
    T, N = len(bubble), len(bubble[0])
    right_phi_x = [np.sum([1 for x in slice if mask > x >= limit]) for slice in bubble]
    time_list = [t if (right_phi_x[t] <= size) else 0 for t in range(T)]
    return next((t for t in time_list[::-1] if t != 0), T-1)

def center_bubble(bubble):
    bubble0 = bubble[0]
    #truncate bubble where it expands relativistically
    limit = phi_upper_bound
    tdecap = time_at_fraction(bubble0, 0.9, np.floor(limit))
    tmax = time_at_fraction(bubble0, 0.95, 4.5)
#    tmax = time_at_fr(bubble0, 0.01, limit)
    T, N = len(bubble0), len(bubble0[0])

    #rotate bubble such that the relativistic bit is perfectly centered
    fld = smoothen(bubble0[tdecap], 1, True)
    kafslj = next(iter([x for x in range(N-1) if fld[x] > limit]))

    if kafslj < 30:
        fld = smoothen(bubble0[tdecap], dx, True)
        vals = [x if (fld[x] > limit and fld[x+1] > limit) else 0 for x in range(N-1)]
        bubbles = [list(g) for k, g in groupby(vals, lambda x: x != 0) if k]
        maxi, imax, keep = 0, 0, 0
        for i in range(len(bubbles)-1):
            maxi = bubbles[i+1][0]-bubbles[i][-1]
            if maxi > imax:
                imax = maxi
                keep = i
    #    first_zero = next((i for i, x in enumerate(vals) if x == 0.), 0) # where bubble begins
        first_zero = int(np.mean([bubbles[keep][0], bubbles[keep][-1]]))
        fld = np.roll(bubble0[tdecap], -first_zero)
        target_peak = int(round(np.mean([x for x in range(N-1) if (fld[x] > limit and fld[x+1] > limit)]))) # center of bubble
    #    print(first_zero, target_peak)    
        angle = int(N//2) - (target_peak + first_zero) # rotation angle needed
        bubble = np.asarray([[np.roll(bubble[col][t], angle) for t in range(tmax)] for col in range(len(bubble))])
    else:
        bubble = np.asarray([bubble[col][:tmax] for col in range(len(bubble))])

    # check which sides the COM travels
    bubble0 = bubble[0]
    tcheck = time_at_fraction(bubble0, 0.02, 2*phi_initial)
#    print(tcheck)
    fld = [smoothen(bubble0[t], dx, True) for t in range(tcheck-100, tcheck, 2)]
#    tv = [np.nanmean([x for x in range(N-1) if (fld[t][x] > 2*phi_initial and fld[t][x+1] > 2*phi_initial)]) for t in range(len(fld))]
    tv = [np.nanmean([x for x in range(N) if fld[t][x] > 2*phi_initial]) for t in range(len(fld))]
    decr = np.sum([1 for i, j in zip(tv, tv[1::]) if i > j])
    incr = np.sum([1 for i, j in zip(tv, tv[1::]) if i < j])
#    print(incr, decr)
    return bubble, decr - incr #np.nanmean(tv) - N//2

def multiply_bubble(bubble, dir, fold):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    bubble = [np.tile(bubble[col], fold) for col in range(len(bubble))] # multiplies bubbles so tail is kept without pbc
    TT, NN = len(bubble[0]), len(bubble[0][0])
    for t in range(TT):
        a, b = int((TT-t)*light_cone) + N - 200, int((TT-t)*light_cone/2) - N//4
        x1, x2 = np.arange(a, NN), np.arange(b)
        if dir < 0:
            x1, x2 = x1 - a, x2 - (b-NN)
        for x in np.append(x1, x2):
            if 0 <= x < NN:
                bubble[0][t][x] = phi_initial
    return np.asarray(bubble)
In [7]:
#all_data1, sims_to_keep1 = np.load(bubbles_file(minSim, 120))
#all_data2, sims_to_keep2 = np.load(bubbles_file(120, 240))
#all_data3, sims_to_keep3 = np.load(bubbles_file(240, 340))
#all_data4, sims_to_keep4 = np.load(bubbles_file(340, nSims))
#all_data = np.concatenate((all_data1, all_data2, all_data3, all_data4)).tolist()
#sims_to_keep = np.concatenate((sims_to_keep1, sims_to_keep2, sims_to_keep3, sims_to_keep4)).tolist()
#del all_data1, all_data2, all_data3, all_data4, sims_to_keep1, sims_to_keep2, sims_to_keep3, sims_to_keep4

# if trying to save the kernel:
all_data, sims_to_keep = np.load(bubbles_file(340, nSims))
In [8]:
normal = [phi_initial, np.mean([sim[1][0] for sim in all_data]), np.mean([sim[2][0] for sim in all_data]), V(phi_initial), 0.5*np.mean([sim[1][0] for sim in all_data])**2]
print(normal)
[3.141592653589793, -6.551177780015445e-05, 1.623650252025356, 0.008, 2.1458965152684045e-09]
In [9]:
#bubble_data = []
#for sim in range(len(all_data))[1:1:]:
#    print(sim)
#    col = 0
#    simulation = all_data[sim]
#    centered, dir = center_bubble(simulation)
#    mult = multiply_bubble(centered, dir, 2)
#    fig, ax = plt.subplots(1, 3, figsize = (12, 3))
#    im0 = ax[0].imshow(simulation[col], aspect='auto', interpolation='none', origin='lower')
#    im1 = ax[1].imshow(centered[col], aspect='auto', interpolation='none', origin='lower')
#    im2 = ax[2].imshow(mult[col], aspect='auto', interpolation='none', origin='lower')
#    clb = plt.colorbar(im0, ax = ax[0]); clb = plt.colorbar(im1, ax = ax[1]); clb = plt.colorbar(im2, ax = ax[2]); plt.show()
#    bubble_data.append(mult)
#np.save(centred_bubbles_file(minSim, nSims), bubble_data)
#check_multiplication = np.load(centred_bubbles_file(minSim, nSims))
#for pic in check_multiplication:
#    plt.figure(0); plt.imshow(pic, aspect='auto', interpolation='none', origin='lower'); plt.show()
In [10]:
def rapidity(vel):
    return np.arctanh(vel)
def gamma(v):
    return ( 1. - v**2. )**(-0.5)

def get_bubble_radius(bubble_slice, filter):
    bubble_slice = smoothen(bubble_slice, filter, True)
    results_half = signal.peak_widths(bubble_slice, signal.find_peaks(bubble_slice)[0], rel_height=0.5)
    return max(results_half[0])/2.

def tanh_pos(x, r0L, r0R, dr, vL, vR):
    wL, wR = dr/gamma(vL), dr/gamma(vR)
    return ( np.tanh( (x - r0L)/wL ) + np.tanh( (r0R - x)/wR ) ) * np.pi/2 + np.pi
def tanh_neg(x, r0L, r0R, dr, vL, vR):
    wL, wR = dr/gamma(vL), dr/gamma(vR)
    return - ( np.tanh( (x - r0L)/wL ) + np.tanh( (r0R - x)/wR ) ) * np.pi/2 + np.pi

def tanh_fit(bubble_slice, axis, prior):
#    plt.plot(axis, bubble_slice, 'ro', axis, [tanh(r, *best_tanh) for r in axis], 'g'); plt.show()
    bounds = ((axis[0], 0, 0, 0, 0), (0, axis[-1], axis[-1], 1, 1))
    if prior is not None:
        if np.mean(bubble_slice) > phi_initial:
            return sco.curve_fit(tanh_pos, axis, bubble_slice, p0=prior, bounds=bounds)[0]
        else:
            return sco.curve_fit(tanh_neg, axis, bubble_slice, p0=prior, bounds=bounds)[0]            
    else:
        if np.mean(bubble_slice) > phi_initial:
            return sco.curve_fit(tanh_pos, axis, bubble_slice, bounds=bounds)[0]
        else:
            return sco.curve_fit(tanh_neg, axis, bubble_slice, bounds=bounds)[0]

def hyperbola1(t, a, b, c):
    return np.sqrt(c + (t - b)**2) + a
def hyperbola2(t, d, e, f):
    return - np.sqrt(f + (t - e)**2) + d

def hypfit_test(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
    try:
        if rr[0] <= rr[-1]:
            fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
            fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
        else:
            fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
            fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
        return True
    except (RuntimeError, TypeError):
        return False
    
def hypfit(tt, rr, qq): # no parameters in common between left and right wall for hyperbolic curve fit
    if rr[0] <= rr[-1]:
        fit1, pcov1 = sco.curve_fit(hyperbola1, tt, rr, p0 = (min(rr), tt[rr.tolist().index(min(rr))], 1e4))
        fit2, pcov2 = sco.curve_fit(hyperbola2, tt, qq, p0 = (max(qq), tt[qq.tolist().index(max(qq))], 1e4))
        traj1, traj2 = hyperbola1(tt, *fit1), hyperbola2(tt, *fit2)
    else:
        fit1, pcov1 = sco.curve_fit(hyperbola2, tt, rr, p0 = (max(rr), tt[rr.tolist().index(max(rr))], 1e4))
        fit2, pcov2 = sco.curve_fit(hyperbola1, tt, qq, p0 = (min(qq), tt[qq.tolist().index(min(qq))], 1e4))
        traj1, traj2 = hyperbola2(tt, *fit1), hyperbola1(tt, *fit2)
#    print('hyperbolic trajectories fit: ', fit1, fit2)
#    plt.plot(rr, tt, 'g-', qq, tt, 'r-', traj1, tt, 'y:', traj2, tt, 'b:') # plot the equation using the fitted parameters
    return tt, traj1, traj2

def get_velocities(tt, rrwallfit, llwallfit):
    #uu = wall travelling along with the com; vv = wall travelling against the com; aa = centre of mass velocity; bb = instantaneous velocity of wall
    uu = intp.splev(tt, intp.splrep(tt, rrwallfit), der=1)
    vv = intp.splev(tt, intp.splrep(tt, llwallfit), der=1)
    for i in range(len(uu)):
        if np.abs(uu[i]) > 1:
            uu[i] = np.sign(uu[i])*(1-np.abs(np.abs(uu[i])-1))
        if np.abs(vv[i]) > 1:
            vv[i] = np.sign(vv[i])*(1-np.abs(np.abs(vv[i])-1))
#    print(uu)
#    print(vv)
    for i in range(len(uu)):
        if str(uu[i]) == 'nan' :
            uu[i] = - np.sign(vv[-1])*0.999
        if str(vv[i]) == 'nan' :
            vv[i] = - np.sign(uu[-1])*0.999
#    print(uu)
#    print(vv)
#    if all(np.abs(i) > 1 for i in uu):
#        uu = np.asarray([np.sign(i)*(1-np.abs(np.abs(i)-1)) for i in uu])
#    if all(np.abs(i) > 1 for i in vv):
#        vv = np.asarray([np.sign(i)*(1-np.abs(np.abs(i)-1)) for i in vv])

    aa = ( 1 + uu*vv - np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / ( uu + vv)
    bb = (-1 + uu*vv + np.sqrt( (-1 + uu**2)*(-1 + vv**2))) / (-uu + vv)
    return uu, vv, aa, bb#, ch1, ch2

def velocity(bubble, bool):
    T, N = len(bubble), len(bubble[0])
    limit = phi_upper_bound - 0.5 #phi_upper_bound# - 0.5*np.abs(phi_upper_bound-np.floor(phi_upper_bound))
    prior = None
    data_list = []
    err = 0

    # find where fit can begin
    window = int(N//10)
    if window > 100: window = 100

    tf = T-1
    endSlice = [i for i in range(N) if np.cos(bubble[tf][i]) > np.cos(limit)]
    if len(endSlice) == 0:
        return 'next'
    while endSlice[0] - window < 0 or endSlice[-1] + window >= N:
        tf -= 2
        endSlice = [i for i in range(N) if np.cos(bubble[tf][i]) > np.cos(limit)]
        if len(endSlice) == 0:
            return 'next'

    for t in range(tf, 0, -1):
        if len(data_list) == 0:
            peaks = endSlice
            target = int(np.round(np.nanmean(peaks)))
            coord_list = np.arange(peaks[0] - window, peaks[-1] + window)
        else:
            xrange = np.arange(data_list[-1][0]-window, data_list[-1][0]+window+1)
            if all(0 <= i < N for i in xrange) and any(np.cos(bubble[t][i]) > np.cos(limit) for i in xrange):
                peaks = [i for i in xrange if np.cos(bubble[t][i]) > np.cos(limit)]
            else:
                break
            target = int(np.round(np.nanmean(peaks)))
            coord_list = np.arange(target - int(np.round(np.abs(data_list[-1][1]))) - window, target + int(np.round(np.abs(data_list[-1][2]))) + window)

        coords = np.asarray([(bubble[t][i%N], int(i-target)) for i in coord_list])
        if err < 5:
            try:
                r0L, r0R, dr, vL, vR = tanh_fit(coords[:,0], coords[:,1], prior)
                prior = r0L, r0R, dr, vL, vR
                data_list.append([target, r0L, r0R, int(t)])
            except (RuntimeError, ValueError, TypeError):
                prior = None
                data_list.append([target, r0L + np.sign(data_list[-1][1]-data_list[0][1])*light_cone*(T/N), r0R + np.sign(data_list[-1][2]-data_list[0][2])*light_cone*(T/N), int(t)])
                err += 1
        else:
            break
    data_list = np.asarray(data_list[::-1])
    targets, r0Ls, r0Rs, time_list = data_list[:,0], data_list[:,1], data_list[:,2], data_list[:,-1]
#    print('tmin, tmax', time_list[0], time_list[-1])

    # get direction of bubble
    radius_diff = np.mean([np.abs(data_list[i,1]) - np.abs(data_list[i,2]) for i in range(len(data_list)//4)])
    if radius_diff > 0:  # if average difference in radius is positive then left radius is on average higher than right radius i.e. bubble travels to the right
        rr, ll = targets + r0Rs, targets + r0Ls
    else:
        rr, ll = targets + r0Ls, targets + r0Rs

    if bool:
        fig, ax0 = plt.subplots(1, 1, figsize = (5, 4))
        ax0.plot(rr, time_list, 'b', ll, time_list, 'y', linewidth='3')
        ax0.set_xlabel('x'); ax0.set_ylabel('t')

    # get velocities from derivative of best fit to wall trajectory
    trunc = 0
    # save copies
    time_list_copy, rr_copy, ll_copy = time_list, rr, ll
    # as lons as it needs..
    while(True):
        # try to fit walls
        if hypfit_test(time_list_copy, rr_copy, ll_copy):
            time_list, rrwallfit, llwallfit = hypfit(time_list_copy, rr_copy, ll_copy)
            uu,vv,aa,bb = get_velocities(time_list, rrwallfit, llwallfit)
            break
        # if fit fails, shorten walls
        else:
            trunc += 25
            if len(time_list_copy) > 100 + trunc*2:
                time_list_copy, rr_copy, ll_copy = time_list_copy[trunc:-trunc:], rr_copy[trunc:-trunc:], ll_copy[trunc:-trunc:]
            if trunc > 200 or len(time_list_copy) < 100:
                return 'next'
            continue
#    if hypfit_test(time_list, rr, ll):
#        time_list, rrwallfit, llwallfit = hypfit(time_list, rr, ll)
#    else:
#        return 'next'
#    uu,vv,aa,bb = get_velocities(time_list, rrwallfit, llwallfit)

#    if bool:
#        fig, [ax0, ax1] = plt.subplots(1, 2, figsize = (15, 4))
#        ax0.plot(rr[-len(time_list):], time_list, 'b', ll[-len(time_list):], time_list, 'y', rrwallfit, time_list, 'r:', llwallfit, time_list, 'r:', linewidth='3')
#        ax0.set_xlabel('x'); ax0.set_ylabel('t')
#        ax1.plot(time_list, uu, 'b:', time_list, vv, 'y:')
#        ax1.plot(time_list, aa, 'r', label=r'v COM')
#        ax1.plot(time_list, bb, 'g', label=r'v walls')
#        ax1.axhline(0, color='darkgray', ls=':')
#        ax1.set_xlabel('t'); ax1.set_ylabel('v(t)'); ax1.legend(); plt.show()
#        plt.savefig(plots_file+'velocity_profile'+str(random.randrange(100))+suffix+'.png');

    list = np.abs(aa - bb)
    if any(str(i) != 'nan' for i in list):
        return -aa[list.tolist().index(np.nanmin(list))]
    else:
        return 'next'
In [11]:
def fold(beta):
    if np.abs(beta) > 0.8:
        return 4
    else:
        return 3

def add_velocities(v1,v2):
    return (v1 + v2) / (1. + v1*v2)

def cut_out_bubble(bubble):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    ext = int(N//8)
    if ext < 100: ext = 200
    xmin = next((i for i in range(N) if bubble0[-1][i] > phi_upper_bound), 0) - ext
    xmax = next((i for i in range(N)[::-1] if bubble0[-1][i] > phi_upper_bound), N-1) + ext
    if xmin < 0: xmin = 0
    if xmax >= N: xmax = N-1
    t, tmin = 0, 0
    while not any(np.cos(bubble0[t][i]) > np.cos(3.17) for i in range(xmin, xmax+1)):
        tmin = t
        t += 1
    return [[[bubble[col][t][x] for x in range(xmin, xmax+1)] for t in range(tmin, T)] for col in range(len(bubble))]

def center_bubble_again(bubble, kind):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    ext = int(N//8)
    if ext < 100: ext = 200
    if kind < 0:
        x1 = 0
        slice = [bubble0[i][x1] for i in range(T)]
        peaks1 = [i for i in slice if i < mask]
        peaks2 = [i for i in range(T) if slice[i] < mask]
        while np.cos(np.nanmean(peaks1)) < np.cos(3.19) and x1 < N-1:
            slice = [bubble0[i][x1] for i in range(T)]
            peaks1 = [i for i in slice if i < mask]
            peaks2 = [i for i in range(T) if slice[i] < mask]
            x1 += 1
        if len(peaks2) > 0:    
            tmax = peaks2[-1]
        else:
            return 'next'
#        print('tmax, x1', tmax, x1)
        xmin = x1 - ext
        if xmin < 0: xmin = 0
        xmax = next((i for i in range(N-1, xmin, -1) if (mask > bubble0[tmax][i] and np.cos(bubble0[tmax][i]) > np.cos(phi_initial))), N-1)
    else:
        x1 = N-1
        slice = [bubble0[i][x1] for i in range(T)]
        peaks1 = [i for i in slice if i < mask]
        peaks2 = [i for i in range(T) if slice[i] < mask]
        while np.cos(np.nanmean(peaks1)) < np.cos(3.19) and x1 > 0:
            slice = [bubble0[i][x1] for i in range(T)]
            peaks1 = [i for i in slice if i < mask]
            peaks2 = [i for i in range(T) if slice[i] < mask]
            x1 -= 1
        if len(peaks2) > 0:    
            tmax = peaks2[-1]
        else:
            return 'next'
#        print('tmax, x1', tmax, x1)
        xmax = x1 + ext
        if xmax > N-1: xmax = N-1
        xmin = next((i for i in range(xmax) if (mask > bubble0[tmax][i] and np.cos(bubble0[tmax][i]) > np.cos(phi_initial))), 0)

    t, tmin = 0, 0
    while not any([(mask > bubble0[t][i] and np.cos(bubble0[t][i]) > np.cos(phi_initial)) for i in range(xmin, xmax+1)]):
        tmin = t
        t += 1

    return [[[bubble[col][t][x] if bubble[col][t][x] != mask else normal[col] for x in range(xmin, xmax+1)] for t in range(tmin, tmax+1)] for col in range(len(bubble))]

def boost(vCOM, ga, x, t):
    return ga * (x - vCOM*t)

def interpolate_bubble(bubble, res):
    NT, N = len(bubble), len(bubble[0])
    bubble = np.asarray(bubble)
    t, x = np.arange(NT), np.arange(N)
    tt, xx = np.meshgrid(t, x)
    f = intp.interp2d(t, x, bubble[tt, xx], kind = 'quintic')
    T, X = np.arange(0, NT, 1/res), np.arange(0, N, 1/res)
    return f(T, X).T

def boost_bubble(bubble, vCOM, res):
    bubble0 = bubble[0]
    T, N = len(bubble0), len(bubble0[0])
    limit = np.ceil(phi_upper_bound)
    ga = gamma(vCOM)

    t0 = int(time_at_size(bubble0, 100, limit)) # arbitrary choice of bubble centre
    x0 = np.mean([i for i in range(N) if mask > bubble0[t0][i] > limit])
    print(t0, x0)
    T0, X0 = boost(vCOM, ga, t0, x0), boost(vCOM, ga, x0, t0)
    deltaT, deltaX = np.abs(T0-t0), np.abs(X0-x0)

    bubble = [interpolate_bubble(bubble[col], res) for col in range(len(bubble))]
    new_bubble = [[[mask for x in range(N)] for t in range(T)] for col in range(len(bubble))]
    for t in range(T):
        for x in range(N):
            tt = int((boost(vCOM, ga, t, x) + np.sign(T//2-t0) * deltaT)*res)
            xx = int((boost(vCOM, ga, x, t) + np.sign(N//2-x0) * deltaX)*res)
            if (T*res > tt >= 0 and N*res > xx >= 0):
                for col in range(len(bubble)):
                    new_bubble[col][t][x] = bubble[col][tt][xx]
    return new_bubble, x0 - N//2
In [12]:
def the_whole_shebang():
    resolution = 3
    for sim in range(len(all_data))[10::]:
        exists = os.path.exists(bubbles_at_rest_file(sim))
        if not exists:

            print(sim)
            bubble, dir = center_bubble(all_data[sim])
            beta = velocity(bubble[0], True)
            if beta == 'next':
                print('sim '+str(sim)+' skipped')
                continue
            elif np.abs(beta) > 0.9:
                beta = np.sign(beta)*0.9
            bubble = multiply_bubble(bubble, dir, fold(beta))

            fig = plt.figure(figsize=(5,5))
            plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
            plt.show()

            res, vtotal = 0, 0
            bool1, bool2 = True, True
            while np.abs(beta) > 0.05:
                bool2 = True
                bubble, dir = boost_bubble(bubble, beta, resolution)
                bubble = center_bubble_again(bubble, dir)
                if bubble == 'next':
                    print('sim '+str(sim)+' interrupted')
                    bool1 = False
                    break

                fig = plt.figure(figsize=(5,5))
                plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
                plt.show()

                prebeta = beta
                vtotal = add_velocities(vtotal, beta)
                print('step done ', beta)

                min_bubble = cut_out_bubble(bubble)
                fig = plt.figure(figsize=(5,5))
                plt.imshow(bubble[0], aspect='auto', interpolation='none', origin='lower')
                plt.show()

                beta = velocity(min_bubble[0], True)
                if beta == 'next':
                    print('sim '+str(sim)+' interrupted')
                    if np.abs(prebeta) >= 0.8:
                        bool1 = False
                    break
                elif np.abs(beta) > 0.9:
                    beta = np.sign(beta)*0.9

#                res += 1

            if not bool2:
                if np.abs(beta) > 0:
                    min_bubble = cut_out_bubble(bubble)
                    vtotal = beta
                    bool2 = True

            if bool1 and bool2:
                print('sim', sim, 'total vel ', vtotal)
                np.save(bubbles_at_rest_file(sim), [min_bubble, vtotal])
    return
In [13]:
#the_whole_shebang()

Average bubble

In [14]:
 def cut_out_coordinates(bubble0):
    T, N = len(bubble0), len(bubble0[0])
    limit = np.floor(phi_upper_bound)
    size = 150 # do not change

    tmin = T-1
    for t in range(T):
        slice = smoothen(bubble0[t], dx, False)
        right_phi = [x for x in range(N-1) if (slice[x] >= limit and slice[x+1] >= limit)]
        if len(right_phi) >= size:
            tmin = t
            break
    if tmin == T-1:
        return None
    if stat.stdev(right_phi) > 100:
        return None
    xcen = int(np.round(np.nanmean(right_phi)))
    xcen_list = [np.round(np.nanmean([x for x in range(N) if bubble0[t][x] >= limit])) for t in range(tmin-50,tmin,2)]
    if stat.stdev(xcen_list) > 50:
        return None

#    fig = plt.figure(figsize=(4,4))
#    plt.plot(xcen, tmin, 'ro')
#    plt.imshow(bubble, aspect='auto', interpolation='none', origin='lower')
#    plt.show()
    ex1 = np.arange(tmin)[::-1]
    ex2 = np.arange(tmin, T)
    ex3 = np.arange(xcen)[::-1]
    ex4 = np.arange(xcen, N)

    if len(ex1) > 500: ex1 = ex1[:500]
    if len(ex2) > 500: ex2 = ex2[:500]
    if len(ex3) > 500: ex3 = ex3[:500]
    if len(ex4) > 500: ex4 = ex4[:500]    
    return ex1, ex2, ex3, ex4

def stack_bubbles(bubble_data):
    extent_list, bad_bubbles = [], []

    for sim in range(len(bubble_data)):
        coords = cut_out_coordinates(bubble_data[sim][0])
        if coords is not None:
            extent_list.append(coords)
        else:
            bad_bubbles.append(sim)

    extent_list = np.asarray(extent_list)
    print('bad bubbles', bad_bubbles)

    bubble_data = [bubble_data[i] for i in range(len(bubble_data)) if i not in bad_bubbles]
    nCols = 5
    for sim in range(len(bubble_data)):
        bubble_data[sim].append([0.5*np.asarray(i)**2 for i in bubble_data[sim][1]])

    uper_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
    upper_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,0][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
    lower_right_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,2][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]
    lower_left_bubble_list = [[[[bubble_data[sim][col][t][x] for t in extent_list[:,1][sim]] for x in extent_list[:,3][sim]] for sim in range(len(bubble_data))] for col in range(nCols)]

    av_mat, av_err_mat = [], []
    for col in range(nCols):
        average_matrix, average_error_matrix = [[],[],[],[]], [[],[],[],[]]
        coord_lists = [uper_right_bubble_list[col], upper_left_bubble_list[col], lower_right_bubble_list[col], lower_left_bubble_list[col]]

        for i in range(4):
            max_lines, max_cols = 500, 500
            for line in range(max_lines):
                average_matrix[i].append([normal[col] for n in range(max_cols)])
                average_error_matrix[i].append([0. for n in range(max_cols)])

            for num_line in range(max_lines):
                for num_col in range(max_cols):
                    meas = []
                    for simulation in coord_lists[i]:
                        if len(simulation) > num_line and len(simulation[0]) > num_col:
                            meas.append(simulation[num_line][num_col])

                    average_matrix[i][num_line][num_col] = np.mean(meas)
                    if len(meas) != 0.:
                        average_error_matrix[i][num_line][num_col] = stat.stdev(meas)/len(meas)
        
        av_mat.append(average_matrix)
        av_err_mat.append(average_error_matrix)
    return [av_mat, av_err_mat], bad_bubbles

def restore_average_bubble(stacked_bubble_data):
    whole_bubbles = []
    for av_bub in stacked_bubble_data:
        whole_bubbles.append([])
        for col in range(len(av_bub)):
            top = np.concatenate((np.flip(np.asarray(np.flip(av_bub[col][0],1)).transpose(),1), np.flip(np.asarray(av_bub[col][1]).transpose(),0)), axis=1)
            bottom = np.concatenate((np.flip(np.asarray(av_bub[col][2]).transpose(),1), np.flip(np.asarray(np.flip(av_bub[col][3],0)).transpose(),1)), axis=1)
            whole_bubbles[-1].append(np.concatenate((top, bottom), axis=0))
    return np.asarray(whole_bubbles)
In [15]:
bubble_list = []
for sim in range(200):#len(all_data)):
    bool = os.path.exists(bubbles_at_rest_file(sim))
    if bool:
        bubble_list.append(np.load(bubbles_at_rest_file(sim))) #deboosted_bubbles2 #deboosted_bubbles_loop_all_sims7
bubble_list = np.asarray(bubble_list)
bubble_list_save = bubble_list[:,0]
bubble_velocities = bubble_list[:,1]
In [16]:
#if testing
#normal = [phi_initial, 0, 0, V(phi_initial), 0]

#either compute average
stacked_bubble_data, bad_bubbles = stack_bubbles(bubble_list_save)
ab, error_ab = restore_average_bubble(stacked_bubble_data)
np.save(average_bubble(minSim, nSims), [ab, error_ab])
#or load it
#ab, error_ab = np.load(average_bubble(minSim, nSims))

#separate columns
bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
error_bubble, error_mom_bubble, error_ge_bubble, error_pe_bubble, error_ke_bubble = error_ab
bad bubbles [3, 5, 10, 12, 14, 16, 26, 27, 28, 29, 35, 40, 45, 46, 51, 54, 57, 61, 62, 63, 64, 66, 70, 73]
In [17]:
#plot discarded bubbles
for i in range(len(bubble_list_save)):
    if i in bad_bubbles:
        fig, ax = plt.subplots(1, 1, figsize=(7,5))
        im = ax.imshow(bubble_list_save[i][0], aspect='auto', interpolation='none', origin='lower')
        clb = plt.colorbar(im, ax = ax)
        plt.show()
In [18]:
bubble_list_save = [bubble_list_save[i] for i in range(len(bubble_list_save)) if i not in bad_bubbles]
bubble_velocities = [bubble_velocities[i] for i in range(len(bubble_velocities)) if i not in bad_bubbles]
check = [bubble_velocities[i] for i in range(len(bubble_velocities)) if i in bad_bubbles]
print('Total number of bubbles available: ', len(bubble_list_save))
Total number of bubbles available:  51
In [19]:
for i in bubble_list_save:
    fig, ax = plt.subplots(1, len(i), figsize=(25,5))
    for col in range(len(i)):
        im = ax[col].imshow(i[col], aspect='auto', interpolation='none', origin='lower')
        clb = plt.colorbar(im, ax = ax[col])
    plt.show()
In [20]:
#plot average bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble
fig, ax = plt.subplots(1, len(ab), figsize=(25, 5))
for col in range(len(ab)):
    im = ax[col].imshow(ab[col], aspect='auto', interpolation='none', origin='lower')
    clb = plt.colorbar(im, ax = ax[col])
plt.savefig(plots_file+'average_bubble'+suffix+'.png')

#plot distribution of boosts
fig, ax = plt.subplots(1, 3, figsize=(15,5))
ax[0].hist(bubble_velocities, bins=len(bubble_velocities))
ax[1].hist([gamma(i) for i in bubble_velocities], bins=len(bubble_velocities))
ax[2].hist([rapidity(i) for i in bubble_velocities], bins=len(bubble_velocities))

[ax[i].set_ylabel('Count') for i in range(len(ax))]
ax[0].set_xlabel(r'$v$'); ax[1].set_xlabel(r'$\gamma$'); ax[2].set_xlabel(r'$w$');
plt.savefig(plots_file+'boost_distrib'+suffix+'.png');
In [21]:
def coleman_match(bubble, colemanSoln):
    minSum = np.sum(colemanSoln)**2.
    tColeman = 0
    for t in range(len(bubble)):
        slice = smoothen(bubble[t], 2*dx, True)
        a = sum([(slice[x] - colemanSoln[x])**2. for x in range(len(slice))])
        if a < minSum:
            minSum, tColeman = a, t
    return tColeman

def plot_coleman_match(bubble, error_bubble, colemanSoln):
    T, N = len(bubble), len(bubble[0])
    colemanSoln = colemanSoln[len(colemanSoln)//2-N//2:len(colemanSoln)//2+N//2]
    tcoleman = coleman_match(bubble, colemanSoln)
    slice, err_slice = bubble[tcoleman], error_bubble[tcoleman]

    fig = plt.figure(figsize=(15,5))
    plt.errorbar(np.arange(N), slice, yerr = err_slice, color = 'k', marker='o', ms=1, alpha = 0.8, ecolor = 'green', label = 'Average bubble')
    plt.plot(np.arange(N), colemanSoln, 'r-', label = 'Best match t = '+str(tcoleman)); plt.legend()
    plt.savefig(plots_file+'average_bubble_profile'+suffix+'.png')
    return tcoleman

def bubble_total_energy(KED, GED, PED):
    T, N = len(KED), len(KED[0])
    KE = [sum(KED[t]) for t in range(T)]
    GE = [sum(GED[t]) for t in range(T)]
    PE = [sum(PED[t]) for t in range(T)]
    TE = [sum(x) for x in zip(KE, GE, PE)]
    return KE, GE, PE, TE

def plot_total_energy_density(KED, GED, PED, tcoleman):
    energies = [KED, GED, PED, KED + GED + PED]
    titles = [r'$\frac{1}{2}(\partial_t \phi)^2$', r'$\frac{1}{2}|\nabla \phi|^2$', r'$V(\phi)$', r'$\frac{1}{2}(\partial_t \phi)^2 +\frac{1}{2}|\nabla \phi|^2+V(\phi)$']

    fig, ax = plt.subplots(2, 2, figsize=(15,10))
    N = len(KED[0])
    for i in range(4):
        if i in [0,1]: j = 0
        else: j = 1
        im = ax[i%2][j].imshow(energies[i], aspect='auto', interpolation='none', origin='lower') 
        clb = plt.colorbar(im, ax = ax[i%2][j]); ax[i%2][j].set_title(titles[i]); ax[i%2][j].set(ylabel=r'$t$')
        ax[i%2][j].plot(np.arange(N), np.ones(N)*tcoleman, 'r-')
    plt.savefig(plots_file+'energy_densities_from_average'+suffix+'.png');
    return

def plot_total_energy(KED, GED, PED, tcoleman):
    energies = bubble_total_energy(KED, GED, PED)
    titles = ['KE', 'GE', 'PE', 'TE']
    ls = ['b-', 'g-', 'y-', 'r-']
    N = len(energies[0])
    fig = plt.figure(figsize=(8,5))
    for i in range(len(energies)):
        plt.plot(np.arange(N), energies[i], ls[i], label = titles[i])
    plt.axhline(np.mean(energies[-1]), color = 'darkgray', ls = '--')
    plt.axvline(tcoleman, color = 'darkgray', ls = '--')
    plt.xlabel('$t$'); plt.legend(); plt.savefig(plots_file+'energies'+suffix+'.png'); plt.show()
    return

def plot_energy_dens_at_t(KED, GED, PED, timeslice):
    TED = [sum(x) for x in zip(KED, GED, PED)]
    energies = [KED, GED, PED, TED]
    titles = ['KED', 'GED', 'PED', 'TED']
    ls = ['b-', 'g-', 'y-', 'r-']
    N = len(KED[0])
    
    fig = plt.figure(figsize=(8,8)); plt.subplots_adjust(hspace=.0)
    gs = gridspec.GridSpec(len(energies), 1, height_ratios=[1]*len(energies))
    ax = [[]]*len(energies); ax[0] = plt.subplot(gs[0])
    for i in range(1, len(ax)): ax[i] = plt.subplot(gs[i], sharex = ax[i-1])
    for i in range(len(ax)):
        ax[i].plot(np.arange(N), energies[i][timeslice], ls[i], label = titles[i])
        if i == 0:
            ax[i].axhline(np.pi, color = 'darkgray', ls='--')
        else:
            ax[i].axhline(np.mean(energies[i][timeslice]), color = 'darkgray', ls='--')            
        ax[i].legend()
    ax[-1].set(xlabel='$x$')
    plt.savefig(plots_file+'energy_densities_at_t'+str(timeslice)+suffix+'.png'); plt.show()
    return
In [22]:
TColeman = plot_coleman_match(bubble, error_bubble, coleman_profile)
plot_total_energy_density(ke_bubble, ge_bubble, pe_bubble, TColeman)
plot_total_energy(ke_bubble, ge_bubble, pe_bubble, TColeman)
#plot_energy_dens_at_t(ke_bubble, ge_bubble, pe_bubble, TColeman)
In [23]:
def GradEnDen(field):
    Nx = len(field)
    NNx = Nx//2+1
    fft_field = scp.fft.fft(field)
    for i in range(len(fft_field)):
        fft_field[i] = 1j*i*dk * fft_field[i]# / Nx
    real_field = scp.fft.ifft(fft_field)
    return [0.5*i.real**2 for i in real_field]

def plot_derived_total_energy_density(field, momentum, filter, t_coleman):
    smooth_field = [smoothen(slice, filter, True) for slice in field]
    smooth_momentum = [smoothen(slice, filter, True) for slice in momentum]
    KED = [0.5*slice**2. for slice in smooth_momentum]
    GED = [GradEnDen(slice) for slice in smooth_field]
    PED = [V(slice) for slice in smooth_field]
    TED = [[KED[t][x]+GED[t][x]+PED[t][x] for x in range(len(KED[0]))] for t in range(len(KED))]
    energies = [KED, GED, PED, TED]
    titles = [r'$\frac{1}{2}(\partial_t \phi)^2$', r'$\frac{1}{2}|\nabla \phi|^2$', r'$V(\phi)$', r'$\frac{1}{2}(\partial_t \phi)^2 +\frac{1}{2}|\nabla \phi|^2+V(\phi)$']

    fig, ax = plt.subplots(2, 2, figsize=(15,10))
    N = len(smooth_field[0])
    for i in range(4):
        if i in [0,1]: j = 0
        else: j = 1
        im = ax[i%2][j].imshow(energies[i], aspect='auto', interpolation='none', origin='lower') 
        clb = plt.colorbar(im, ax = ax[i%2][j]); ax[i%2][j].set_title(titles[i]); ax[i%2][j].set(ylabel=r'$t$')
        ax[i%2][j].plot(np.arange(N), np.ones(N)*t_coleman, 'r-')
    plt.savefig(plots_file+'energy_densities'+suffix+'.png');
    return

def plot_TE_with_Nbubbles(tcoleman):
    bubble_averages_list = []
    for i in range(Nbubs):
        bool = os.path.exists(average_of_N_bubbles(i+1))
        if bool:
            bubble_averages_list.append(np.load(average_of_N_bubbles(i+1)))

    fig = plt.figure(figsize=(8,6))
    for bubble_average in bubble_averages_list:
        ab, error_ab, length = bubble_average
        bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
        energies = bubble_total_energy(ke_bubble, ge_bubble, pe_bubble)
        plt.plot(np.arange(len(energies[-1])), energies[-1], label = length)
        plt.axhline(np.mean(energies[-1]), color = 'darkgray', ls = '--')
    plt.axvline(tcoleman, color = 'darkgray', ls = '--')
    plt.xlabel('$t$'); plt.legend()
    plt.savefig(plots_file+'TE_vs_nbubbles1'+suffix+'.png');

    fig = plt.figure(figsize=(8,6))
    for bubble_average in bubble_averages_list:
        ab, error_ab, length = bubble_average
        bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
        GED, PED, KED = ge_bubble[tcoleman], pe_bubble[tcoleman], ke_bubble[tcoleman]
        TED = GED + PED + KED
        plt.plot(length, sum(TED), 'ro')#, np.mean(TED))
    plt.xlabel('$N$'); plt.ylabel('$E$'); plt.legend()
    plt.savefig(plots_file+'TE_vs_nbubbles2'+suffix+'.png');
    return

def plot_TE_with_filter(tcoleman):
    bubble_average = np.load(average_of_N_bubbles(Nbubs))
    ab, error_ab, length = bubble_average
    bubble, mom_bubble, ge_bubble, pe_bubble, ke_bubble = ab
    slice = bubble[tcoleman]
    mom_slice = mom_bubble[tcoleman]

    fig = plt.figure(figsize=(8,6))
    for filter in np.linspace(dx, 1000*dx, 200):
        smooth_field = smoothen(slice, filter, True)
        smooth_momentum = smoothen(mom_slice, filter, True)
        KED = 0.5*smooth_momentum**2.
        GED = GradEnDen(smooth_field)
        PED = V(smooth_field)
        TED = KED + GED + PED
        plt.plot(filter, np.sum(TED), 'ro', ms=2)#, 'o', label='{:.3f}'.format(filter))
    plt.axhline(0, ls='--', color='darkgray')
    plt.xlabel(r'$\sigma$'); plt.ylabel('$E$'); plt.legend()
    plt.savefig(plots_file+'TE_vs_filter'+suffix+'.png');
    return
In [24]:
plot_derived_total_energy_density(bubble, mom_bubble, 5, TColeman)
In [25]:
#average x bubbles at a time; save
Nbubs = len(bubble_list_save)
step = 4
for i in range(0, Nbubs, step):
    qq = len(bubble_list_save)
    bubble_puzzle, _ = stack_bubbles(bubble_list_save)
    ab, error_ab = restore_average_bubble(bubble_puzzle)
    np.save(average_of_N_bubbles(qq), [ab, error_ab, qq])
    for i in range(step):
        if len(bubble_list_save) > 0:
            del bubble_list_save[0]
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
bad bubbles []
In [26]:
plot_TE_with_Nbubbles(TColeman)
No handles with labels found to put in legend.
In [27]:
#Nbubs = 62
plot_TE_with_filter(TColeman)
No handles with labels found to put in legend.
In [ ]: